Simple CJS model integrated with a growth in weight model to get phi, p, and growth estimates to develop a production model.
In all the models below, 1 = not observed and 2 = observed in the input encounter data.
Encounter data are available here in the eh.csv file. Weight data are in weight.csv
Model code is in ./models/production/modelProduction_OB_functionsToSource.R
Model is run ‘by hand’ in ./models/modelProduction_OB_makeFile.R
Functions for this qmd file in ./models/qmdProduction_OB_functionsToSource.R
There are 3 models here: All cohorts:
model 1, with quadratic weight effect on phi By cohort:
model 2, with linear weight effect on phi model 3, with quadratic weight effect on phi
all <- tar_read(cdWB_electro_target)
table(all %>% filter(river == "wb obear")|> dplyr::select(ageInSamples))ageInSamples
0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15
1308 1890 1007 837 1486 875 446 298 444 159 155 57 47 17 16 3
16 18
3 3 cohorts <- 2002:2014Using model #3 from modelsCMR_Growth_NN_OB.qmd as a staring point for the models, but adapting the model by
Prob not 1) Extending AgeInSamples from 1-11 to 1-x to allow bigger fish to be present for the production estimates.
Prob not 2) Loop over first[i]:lastAIS so fish have the chance to survive to large size.
Model with phi and p for each age-in-samples (t = column in the encounter history file) and individual (i)
Probability of survival (phi) model structure:
logit(phi[i,t]) <-
phiInt[i,t] +
phiBeta[1,i,t] * weight[i,t] Probability of capture (p) model structure:
p(t,i) <- pInt(t,i)
Growth rate (gr) model structure:
gr(t,i) <- grInt(t)
Survival/growth rate interaction:
Multiplicative
Model runs:
library(targets)
library(nimble)
# Following https://oliviergimenez.github.io/bayesian-cr-workshop/worksheets/4_demo.html
#
#
# To get the mMCMCSummaryRMNA funcion which I adapted to deal with NAs
source('./models/production/modelProduction_OB_functionsToSource.R')
source('./models/production/qmdProduction_OB_functionsToSource.R')
modelNum <- 2 # phi * growth# all cohorts have the same code for a given model. show model code for one of them here
load('./models/production/runsOut/production_OB_model_2_2002_mostRecent.RData')
d$modelCode{
for (i in 1:N) {
weight[i, first[i]] ~ dnorm(meanWeight_AISstd[first[i]],
sd = 2)
weightDATAstd[i, first[i]] ~ dnorm(weight[i, first[i]],
sd = weightSD)
for (t in (first[i]):(last[i] - 1)) {
weight[i, t + 1] <- weight[i, t] + gr[i, t] * sampleIntervalDATA[i,
t]/90
weightDATAstd[i, t + 1] ~ dnorm(weight[i, t + 1],
sd = weightSD)
}
}
for (i in 1:N) {
for (t in first[i]:(last[i] - 1)) {
gr[i, t] ~ dnorm(grIntT[t], sd = 0.5)
}
}
for (t in 1:(T - 1)) {
grIntT[t] ~ dnorm(0, sd = 1000)
}
delta[1] <- 1
delta[2] <- 0
for (i in 1:N) {
for (t in first[i]:(last[i] - 1)) {
gamma[1, 1, t, i] <- phi[i, t]
gamma[1, 2, t, i] <- 1 - phi[i, t]
gamma[2, 1, t, i] <- 0
gamma[2, 2, t, i] <- 1
}
}
for (t in 1:(T - 1)) {
p[t] ~ dunif(0, 1)
omega[1, 1, t] <- 1 - p[t]
omega[1, 2, t] <- p[t]
omega[2, 1, t] <- 1
omega[2, 2, t] <- 0
}
for (i in 1:N) {
for (t in first[i]:(last[i] - 1)) {
logit(phi[i, t]) <- phiInt[i, t] + phiBeta[1, i,
t] * weight[i, t]
}
}
for (i in 1:N) {
for (t in first[i]:(last[i] - 1)) {
phiInt[i, t] ~ dnorm(phiIntT[t], sd = 2)
for (b in 1:1) {
phiBeta[b, i, t] ~ dnorm(phiBetaT[b, t], sd = 1/0.667)
}
}
}
for (t in 1:(T - 1)) {
phiIntT[t] ~ T(dnorm(0, sd = 0.667), -3.5, 3.5)
phiIntTOut[t] <- ilogit(phiIntT[t])
phiBetaT[1, t] ~ dnorm(0, sd = 1/0.667)
}
for (i in 1:N) {
for (t in 1:T) {
weightZ01[i, t] <- weightIfAvailable(weight[i, t],
sdWeight, meanWeight, z[i, t], available[i, t])
z01[i, t] <- aliveIfAvailable(z[i, t], available[i,
t])
}
}
for (t in 1:T) {
sumWeightZ01[t] <- sum(weightZ01[1:N, t])
sumZ01[t] <- sum(z01[1:N, t])
}
for (i in 1:N) {
z[i, first[i]] ~ dcat(delta[1:2])
for (j in first[i]:(last[i] - 1)) {
z[i, j + 1] ~ dcat(gamma[z[i, j], 1:2, j, i])
yDATA[i, j + 1] ~ dcat(omega[z[i, j + 1], 1:2, j])
}
}
}getSummaryByCohort <- function(cohort, modelNum = 3) {
# Load the data
load(paste0('./models/production/runsOut/production_OB_model', '_', modelNum, '_', cohort, '_mostRecent.RData'))
# Input data
eh <- tar_read_raw(paste0('eh_OB_', cohort, '_target'))
print(paste0('Run time = ', round(d$runTime, 3), ' ', attr(d$runTime, "units")))
MCMCplot(object = d$mcmc, params = c("sumZ01"), ref_ovl = TRUE)
MCMCplot(object = d$mcmc, params = c("sumWeightZ01"), ref_ovl = TRUE)
MCMCplot(object = d$mcmc, params = c("p"), ref_ovl = TRUE)
MCMCplot(object = d$mcmc, params = c("grIntT"), ref_ovl = TRUE)
MCMCplot(object = d$mcmc, params = c("phiIntTOut"), ref_ovl = TRUE)
MCMCplot(object = d$mcmc, params = c("phiBetaT"), ref_ovl = TRUE)
get_MCMCTrace_OB(d)
# Create the summary list
summary_list <- list(
eh = eh,
stats = kable(as_tibble(d$runData), caption = "Run statistics"),
runTime = paste0('Run time = ', round(d$runTime, 3), ' ', attr(d$runTime, "units")),
#traces = get_MCMCTrace_OB(d),
summaryMod3_tt_growth = MCMCsummary(object = d$mcmc, params = c(
#"phiIntT",
"phiIntTOut",
"p",
"grIntT",
"phiBetaT",
"sumWeightZ01",
"sumZ01"
), round = 3) %>%
rownames_to_column(var = "var"),
summary_OB = get_summary_OB(d, eh)
)
return(summary_list)
}
#s=summary_OB[, order(summary_OB$ind, summary_OB$t)]s2002_2 = getSummaryByCohort(2002, modelNum)[1] "Run time = 21.474 mins"
s2003_2 = getSummaryByCohort(2003, modelNum)[1] "Run time = 27.403 mins"
s2004_2 = getSummaryByCohort(2004, modelNum)[1] "Run time = 17.544 mins"
s2005_2 = getSummaryByCohort(2005, modelNum)[1] "Run time = 8.699 mins"
s2006_2 = getSummaryByCohort(2006, modelNum)[1] "Run time = 14.192 mins"
s2007_2 = getSummaryByCohort(2007, modelNum)[1] "Run time = 11.018 mins"
s2008_2 = getSummaryByCohort(2008, modelNum)[1] "Run time = 9.901 mins"
s2009_2 = getSummaryByCohort(2009, modelNum)[1] "Run time = 55.498 mins"
s2010_2 = getSummaryByCohort(2010, modelNum)[1] "Run time = 7.013 mins"
s2011_2 = getSummaryByCohort(2011, modelNum)[1] "Run time = 6.584 mins"
s2012_2 = getSummaryByCohort(2012, modelNum)[1] "Run time = 1.09 hours"
s2013_2 = getSummaryByCohort(2013, modelNum)[1] "Run time = 8.649 mins"
s2014_2 = getSummaryByCohort(2014, modelNum)[1] "Run time = 6.429 mins"
sAll_2 <- tibble(.rows = 0) |>
add_column(!!!s2002_2$summary_OB[0,])
# Fill tibble
for (cohort in cohorts) {
summary_obj_2 <- get(paste0("s", cohort, "_2"))$summary_OB
summary_obj_2$cohort <- cohort
sAll_2 <- bind_rows(sAll_2, summary_obj_2)
}
sAll_2$cohort_ind <- paste0(sAll_2$cohort, "_", sAll_2$ind)#ojs_define(numTags_OJS_mod3 = dim(s2002$eh$tags)[1]) # all cohorts have the same eh
ojs_define(summary_OB_OJS_all_2 = sAll_2)
#ojs_define(summary_tt_z_OJS = (summary_tt_z_mod3))summary_OB_OJS_all_T_2 = transpose(summary_OB_OJS_all_2)import { range } from "d3";
viewof selected_cohort_2 = Inputs.select(d3.range(2002, 2015), {
label: "Which cohort?",
value: 2002,
step: 1,
multiple: true
})summary_OB_OJS_all_T_selected_cohort_2 = summary_OB_OJS_all_T_2.filter((d) =>
(selected_cohort_2.includes(d.cohort))
)uniqueInds_2 = [...new Set(summary_OB_OJS_all_T_selected_cohort_2.map(d => d.cohort_ind))].sort();
viewof selectInd_mod_2 = Inputs.select(uniqueInds_2, {
label: "Which fish?",
value: 1,
step: 1,
multiple: true
})
summary_OB_OJS_all_T_selected_cohort_ind_2 = summary_OB_OJS_all_T_selected_cohort_2.filter((d) =>
(selectInd_mod_2.includes(d.cohort_ind))
)Black dots/line is estimated mass and orange dots are observed masses. The green line is the first observation and the red line is the last observation. The number in the upper right corner of each panel is the fish’s cohort.
Black dots/line is estimated probability of survival and orange dots are encountered yes (y = 0.8)/no (y = 0). The green line is the first observation and the red line is the last observation. The number in the upper right corner of each panel is the fish’s cohort.
Plot.plot({
width: width,
height: 350,
inset: 10,
color: {
scheme: "greys"
},
x: { label: "Age/season combination" },
y: { label: "Probability of survival" },
marks: [
Plot.frame(),
Plot.dot(summary_OB_OJS_all_T_selected_cohort_ind_2, {
x: "t",
y: "pSurv"
}),
Plot.dot(summary_OB_OJS_all_T_selected_cohort_ind_2, {
x: "t",
y: "enc8",
fill: "orange"
}),
Plot.line(summary_OB_OJS_all_T_selected_cohort_ind_2, {
x: "t",
y: "pSurv"
}),
Plot.ruleX(summary_OB_OJS_all_T_selected_cohort_ind_2, {
x: "first",
y: 1,
"stroke": "green"
}),
Plot.ruleX(summary_OB_OJS_all_T_selected_cohort_ind_2, {
x: "last",
y: 1,
"stroke": "red"
})
,
Plot.text(summary_OB_OJS_all_T_selected_cohort_ind_2,
Plot.selectFirst({
x: 11,
y: 1,
text: d => d.cohort
})
)
],
facet: {
data: summary_OB_OJS_all_T_selected_cohort_ind_2,
x: "ind"
}
})Plot.plot({
width: width,
height: 350,
inset: 10,
color: {
scheme: "greys"
},
x: { label: "Age/season combination" },
y: { label: "Standardized body mass (mg)" },
marks: [
Plot.frame(),
Plot.dot(summary_OB_OJS_all_T_selected_cohort_ind_2, {
x: "t",
y: "mean_weight"
}),
Plot.line(summary_OB_OJS_all_T_selected_cohort_ind_2, {
x: "t",
y: "mean_weight"
}),
<!-- Plot.dot(d, { -->
<!-- x: "t", -->
<!-- y: "mean_gr", -->
<!-- fill: "green" -->
<!-- }), -->
Plot.dot(summary_OB_OJS_all_T_selected_cohort_ind_2, {
x: "t",
y: "weightDATA_std",
fill: "orange"
}),
Plot.ruleX(summary_OB_OJS_all_T_selected_cohort_ind_2, {
x: "first",
y: 3,
"stroke": "green"
}),
Plot.ruleX(summary_OB_OJS_all_T_selected_cohort_ind_2, {
x: "last",
y: 3,
"stroke": "red"
}),
Plot.text(summary_OB_OJS_all_T_selected_cohort_ind_2,
Plot.selectFirst({
x: 1,
y: 4,
frameAnchor: "top-left",
text: (d) => "Cohort = " + d.cohort
})
),
Plot.text(summary_OB_OJS_all_T_selected_cohort_ind_2,
Plot.selectFirst({
x: 1,
y: 3.5,
frameAnchor: "top-left",
text: (d) => "Residual = " + d.meanResid
})
)
],
facet: {
data: summary_OB_OJS_all_T_selected_cohort_ind_2,
x: "ind"
}
})Combine scohort$summaryMod2_tt_growth data into one data frame
sParamsAll0_2 <- tibble(.rows = 0) |>
add_column(!!!s2002_2$summaryMod3_tt_growth[0,]) # need to change to mod2
# Fill tibble
for (cohort in cohorts) {
summary_obj_2 <- get(paste0("s", cohort, "_2"))$summaryMod3_tt_growth # need to change to mod2
summary_obj_2$cohort <- cohort
sParamsAll0_2 <- bind_rows(sParamsAll0_2, summary_obj_2)
}
# remove brackets for filtering
sParamsAll0_2$varNoIndex <- sParamsAll0_2$var |>
str_remove("\\[.*\\]")
# for numeric ordering
sParamsAll1_2 <- sParamsAll0_2 |> filter(varNoIndex != "phiBetaT") |>
mutate(
varNumeric1 = var |>
str_extract("\\[(\\d+)\\]") %>%
str_remove_all("\\[|\\]") %>%
as.numeric()
)
sParamsAll2_2 <- sParamsAll0_2 |> filter(varNoIndex == "phiBetaT") |>
mutate(
varNumeric1 = var |>
str_extract("\\[(\\d+),") %>% # Changed to match number before comma
str_remove_all("\\[|,") %>% # Remove [ and comma
as.numeric(),
varNumeric2 = var |>
str_extract(",\\s*(\\d+)\\]") %>%
str_remove_all(",|\\s|\\]") %>%
as.numeric()
)
sParamsAll_2 = bind_rows(sParamsAll1_2, sParamsAll2_2)ggplot(sParamsAll_2 |> filter(varNoIndex == "p"), aes(varNumeric1, mean, group = cohort, color = cohort)) +
geom_point() +
geom_line() +
xlab("AgeInSamples")
ggplot(sParamsAll_2 |> filter(varNoIndex == "phiIntTOut"), aes(varNumeric1, mean, group = cohort, color = factor(cohort))) +
geom_point() +
geom_line() +
xlab("AgeInSamples")
ggplot(sParamsAll_2 |> filter(varNoIndex == "grIntT"), aes(varNumeric1, mean, group = cohort, color = factor(cohort))) +
geom_point() +
geom_line() +
xlab("AgeInSamples")
ggplot(sParamsAll_2 |> filter(varNoIndex == "grIntT", mean > -5, mean < 5), aes(varNumeric1, mean, group = cohort, color = factor(cohort))) +
geom_point() +
geom_line() +
xlab("AgeInSamples")
ggplot(sParamsAll_2 |> filter(varNoIndex == "phiBetaT", mean > -5, mean < 5), aes(varNumeric2, mean, group = cohort, color = factor(cohort))) +
geom_point() +
geom_line() +
xlab("AgeInSamples") +
facet_wrap((~varNumeric1))
ggplot(sParamsAll_2 |> filter(varNoIndex == "sumWeightZ01"), aes(varNumeric1, mean, group = cohort, color = factor(cohort))) +
geom_point() +
geom_line() +
xlab("AgeInSamples")
ggplot(sParamsAll_2 |> filter(varNoIndex == "sumZ01"), aes(varNumeric1, mean, group = cohort, color = factor(cohort))) +
geom_point() +
geom_line() +
xlab("AgeInSamples")
Need to get rid of squared term for cohort-specific models
phiInt_2 <- sParamsAll_2 |> filter(varNoIndex == "phiIntTOut") |>
dplyr::select(mean, varNumeric1, cohort) |>
rename(phiIntOut = mean, ais = varNumeric1)
phiBetas_2 <- sParamsAll_2 |> filter(varNoIndex == "phiBetaT", mean > -5, mean < 5) |>
dplyr::select(mean, varNumeric1, varNumeric2, cohort) |>
pivot_wider(names_from = varNumeric1, values_from = mean, names_prefix = "beta") |>
rename(ais = varNumeric2)
phiSize_2 = expand.grid(weight = seq(-2,2,0.25), ais = 1:11, cohort = cohorts) |>
left_join(phiInt_2) |>
left_join(phiBetas_2) |>
mutate(pSurv = ilogit(phiIntOut + beta1*weight))# + beta_2*weight^2))
ggplot(phiSize_2, aes(weight, pSurv, color = factor(cohort), group = cohort)) +
geom_line() +
geom_point() +
facet_wrap(~ais)
Using model #3 from modelsCMR_Growth_NN_OB.qmd as a staring point for the models, but adapting the model by
Prob not 1) Extending AgeInSamples from 1-11 to 1-x to allow bigger fish to be present for the production estimates.
Prob not 2) Loop over first[i]:lastAIS so fish have the chance to survive to large size.
Model with phi and p for each age-in-samples (t = column in the encounter history file) and individual (i)
Probability of survival (phi) model structure:
logit(phi[i,t]) <-
phiInt[i,t] +
phiBeta[1,i,t] * weight[i,t] +
phiBeta[2,i,t] * weight[i,t] * weight[i,t] Probability of capture (p) model structure:
p(t,i) <- pInt(t,i)
Growth rate (gr) model structure:
gr(t,i) <- grInt(t)
Survival/growth rate interaction:
Multiplicative
Model runs:
library(targets)
library(nimble)
# Following https://oliviergimenez.github.io/bayesian-cr-workshop/worksheets/4_demo.html
#
#
# To get the mMCMCSummaryRMNA funcion which I adapted to deal with NAs
#source('./models/production/modelProduction_OB_functionsToSource.R')
#source('./models/production/qmdProduction_OB_functionsToSource.R')
modelNum <- 3 # phi * growth# all cohorts have the same code for a given model. show model code for one of them here
load('./models/production/runsOut/production_OB_model_3_2002_mostRecent.RData')
d$modelCode{
for (i in 1:N) {
weight[i, first[i]] ~ dnorm(meanWeight_AISstd[first[i]],
sd = 2)
weightDATAstd[i, first[i]] ~ dnorm(weight[i, first[i]],
sd = weightSD)
for (t in (first[i]):(last[i] - 1)) {
weight[i, t + 1] <- weight[i, t] + gr[i, t] * sampleIntervalDATA[i,
t]/90
weightDATAstd[i, t + 1] ~ dnorm(weight[i, t + 1],
sd = weightSD)
}
}
for (i in 1:N) {
for (t in first[i]:(last[i] - 1)) {
gr[i, t] ~ dnorm(grIntT[t], sd = 0.5)
}
}
for (t in 1:(T - 1)) {
grIntT[t] ~ dnorm(0, sd = 1000)
}
delta[1] <- 1
delta[2] <- 0
for (i in 1:N) {
for (t in first[i]:(last[i] - 1)) {
gamma[1, 1, t, i] <- phi[i, t]
gamma[1, 2, t, i] <- 1 - phi[i, t]
gamma[2, 1, t, i] <- 0
gamma[2, 2, t, i] <- 1
}
}
for (t in 1:(T - 1)) {
p[t] ~ dunif(0, 1)
omega[1, 1, t] <- 1 - p[t]
omega[1, 2, t] <- p[t]
omega[2, 1, t] <- 1
omega[2, 2, t] <- 0
}
for (i in 1:N) {
for (t in first[i]:(last[i] - 1)) {
logit(phi[i, t]) <- phiInt[i, t] + phiBeta[1, i,
t] * weight[i, t] + phiBeta[2, i, t] * weight[i,
t] * weight[i, t]
}
}
for (i in 1:N) {
for (t in first[i]:(last[i] - 1)) {
phiInt[i, t] ~ dnorm(phiIntT[t], sd = 2)
for (b in 1:2) {
phiBeta[b, i, t] ~ dnorm(phiBetaT[b, t], sd = 1/0.667)
}
}
}
for (t in 1:(T - 1)) {
phiIntT[t] ~ T(dnorm(0, sd = 0.667), -3.5, 3.5)
phiIntTOut[t] <- ilogit(phiIntT[t])
phiBetaT[1, t] ~ dnorm(0, sd = 1/0.667)
phiBetaT[2, t] ~ dnorm(0, sd = 1/0.667)
}
for (i in 1:N) {
for (t in 1:T) {
weightZ01[i, t] <- weightIfAvailable(weight[i, t],
sdWeight, meanWeight, z[i, t], available[i, t])
z01[i, t] <- aliveIfAvailable(z[i, t], available[i,
t])
}
}
for (t in 1:T) {
sumWeightZ01[t] <- sum(weightZ01[1:N, t])
sumZ01[t] <- sum(z01[1:N, t])
}
for (i in 1:N) {
z[i, first[i]] ~ dcat(delta[1:2])
for (j in first[i]:(last[i] - 1)) {
z[i, j + 1] ~ dcat(gamma[z[i, j], 1:2, j, i])
yDATA[i, j + 1] ~ dcat(omega[z[i, j + 1], 1:2, j])
}
}
}s2002_3 = getSummaryByCohort(2002, modelNum)[1] "Run time = 32.092 mins"
s2003_3 = getSummaryByCohort(2003, modelNum)[1] "Run time = 37.935 mins"
s2004_3 = getSummaryByCohort(2004, modelNum)[1] "Run time = 28.262 mins"
s2005_3 = getSummaryByCohort(2005, modelNum)[1] "Run time = 12.692 mins"
s2006_3 = getSummaryByCohort(2006, modelNum)[1] "Run time = 20.456 mins"
s2007_3 = getSummaryByCohort(2007, modelNum)[1] "Run time = 15.911 mins"
s2008_3 = getSummaryByCohort(2008, modelNum)[1] "Run time = 14.706 mins"
s2009_3 = getSummaryByCohort(2009, modelNum)[1] "Run time = 1.453 hours"
s2010_3 = getSummaryByCohort(2010, modelNum)[1] "Run time = 10.1 mins"
s2011_3 = getSummaryByCohort(2011, modelNum)[1] "Run time = 9.941 mins"
s2012_3 = getSummaryByCohort(2012, modelNum)[1] "Run time = 1.39 hours"
s2013_3 = getSummaryByCohort(2013, modelNum)[1] "Run time = 13.827 mins"
s2014_3 = getSummaryByCohort(2014, modelNum)[1] "Run time = 8.986 mins"
sAll_3 <- tibble(.rows = 0) |>
add_column(!!!s2002_3$summary_OB[0,])
# Fill tibble
for (cohort in cohorts) {
summary_obj_3 <- get(paste0("s", cohort, "_3"))$summary_OB
summary_obj_3$cohort <- cohort
sAll_3 <- bind_rows(sAll_3, summary_obj_3)
}
sAll_3$cohort_ind <- paste0(sAll_3$cohort, "_", sAll_3$ind)#ojs_define(numTags_OJS_mod3 = dim(s2002$eh$tags)[1]) # all cohorts have the same eh
ojs_define(summary_OB_OJS_all_3 = sAll_3)
#ojs_define(summary_tt_z_OJS = (summary_tt_z_mod3))summary_OB_OJS_all_T_3 = transpose(summary_OB_OJS_all_3)import { range } from "d3";
viewof selected_cohort_3 = Inputs.select(d3.range(2002, 2015), {
label: "Which cohort?",
value: 2002,
step: 1,
multiple: true
})summary_OB_OJS_all_T_selected_cohort_3 = summary_OB_OJS_all_T_3.filter((d) =>
(selected_cohort_3.includes(d.cohort))
)uniqueInds_3 = [...new Set(summary_OB_OJS_all_T_selected_cohort_3.map(d => d.cohort_ind))].sort();
viewof selectInd_mod_3 = Inputs.select(uniqueInds_3, {
label: "Which fish?",
value: 1,
step: 1,
multiple: true
})
summary_OB_OJS_all_T_selected_cohort_ind_3 = summary_OB_OJS_all_T_selected_cohort_3.filter((d) =>
(selectInd_mod_3.includes(d.cohort_ind))
)Black dots/line is estimated mass and orange dots are observed masses. The green line is the first observation and the red line is the last observation. The number in the upper right corner of each panel is the fish’s cohort.
Black dots/line is estimated probability of survival and orange dots are encountered yes (y = 0.8)/no (y = 0). The green line is the first observation and the red line is the last observation. The number in the upper right corner of each panel is the fish’s cohort.
Plot.plot({
width: width,
height: 350,
inset: 10,
color: {
scheme: "greys"
},
x: { label: "Age/season combination" },
y: { label: "Probability of survival" },
marks: [
Plot.frame(),
Plot.dot(summary_OB_OJS_all_T_selected_cohort_ind_3, {
x: "t",
y: "pSurv"
}),
Plot.dot(summary_OB_OJS_all_T_selected_cohort_ind_3, {
x: "t",
y: "enc8",
fill: "orange"
}),
Plot.line(summary_OB_OJS_all_T_selected_cohort_ind_3, {
x: "t",
y: "pSurv"
}),
Plot.ruleX(summary_OB_OJS_all_T_selected_cohort_ind_3, {
x: "first",
y: 1,
"stroke": "green"
}),
Plot.ruleX(summary_OB_OJS_all_T_selected_cohort_ind_3, {
x: "last",
y: 1,
"stroke": "red"
})
,
Plot.text(summary_OB_OJS_all_T_selected_cohort_ind_3,
Plot.selectFirst({
x: 11,
y: 1,
text: d => d.cohort
})
)
],
facet: {
data: summary_OB_OJS_all_T_selected_cohort_ind_3,
x: "ind"
}
})Plot.plot({
width: width,
height: 350,
inset: 10,
color: {
scheme: "greys"
},
x: { label: "Age/season combination" },
y: { label: "Standardized body mass (mg)" },
marks: [
Plot.frame(),
Plot.dot(summary_OB_OJS_all_T_selected_cohort_ind_3, {
x: "t",
y: "mean_weight"
}),
Plot.line(summary_OB_OJS_all_T_selected_cohort_ind_3, {
x: "t",
y: "mean_weight"
}),
<!-- Plot.dot(d, { -->
<!-- x: "t", -->
<!-- y: "mean_gr", -->
<!-- fill: "green" -->
<!-- }), -->
Plot.dot(summary_OB_OJS_all_T_selected_cohort_ind_3, {
x: "t",
y: "weightDATA_std",
fill: "orange"
}),
Plot.ruleX(summary_OB_OJS_all_T_selected_cohort_ind_3, {
x: "first",
y: 3,
"stroke": "green"
}),
Plot.ruleX(summary_OB_OJS_all_T_selected_cohort_ind_3, {
x: "last",
y: 3,
"stroke": "red"
}),
Plot.text(summary_OB_OJS_all_T_selected_cohort_ind_3,
Plot.selectFirst({
x: 1,
y: 4,
frameAnchor: "top-left",
text: (d) => "Cohort = " + d.cohort
})
),
Plot.text(summary_OB_OJS_all_T_selected_cohort_ind_3,
Plot.selectFirst({
x: 1,
y: 3.5,
frameAnchor: "top-left",
text: (d) => "Residual = " + d.meanResid
})
)
],
facet: {
data: summary_OB_OJS_all_T_selected_cohort_ind_3,
x: "ind"
}
})Combine scohort$summaryMod2_tt_growth data into one data frame
sParamsAll0_3 <- tibble(.rows = 0) |>
add_column(!!!s2002_3$summaryMod3_tt_growth[0,])
# Fill tibble
for (cohort in cohorts) {
summary_obj_3 <- get(paste0("s", cohort, "_3"))$summaryMod3_tt_growth
summary_obj_3$cohort <- cohort
sParamsAll0_3 <- bind_rows(sParamsAll0_3, summary_obj_3)
}
# remove brackets for filtering
sParamsAll0_3$varNoIndex <- sParamsAll0_3$var |>
str_remove("\\[.*\\]")
# for numeric ordering
sParamsAll1_3 <- sParamsAll0_3 |> filter(varNoIndex != "phiBetaT") |>
mutate(
varNumeric1 = var |>
str_extract("\\[(\\d+)\\]") %>%
str_remove_all("\\[|\\]") %>%
as.numeric()
)
sParamsAll2_3 <- sParamsAll0_3 |> filter(varNoIndex == "phiBetaT") |>
mutate(
varNumeric1 = var |>
str_extract("\\[(\\d+),") %>% # Changed to match number before comma
str_remove_all("\\[|,") %>% # Remove [ and comma
as.numeric(),
varNumeric2 = var |>
str_extract(",\\s*(\\d+)\\]") %>%
str_remove_all(",|\\s|\\]") %>%
as.numeric()
)
sParamsAll_3 = bind_rows(sParamsAll1_3, sParamsAll2_3)ggplot(sParamsAll_3 |> filter(varNoIndex == "p"), aes(varNumeric1, mean, group = cohort, color = cohort)) +
geom_point() +
geom_line() +
xlab("AgeInSamples")
ggplot(sParamsAll_3 |> filter(varNoIndex == "phiIntTOut"), aes(varNumeric1, mean, group = cohort, color = cohort)) +
geom_point() +
geom_line() +
xlab("AgeInSamples")
ggplot(sParamsAll_3 |> filter(varNoIndex == "grIntT"), aes(varNumeric1, mean, group = cohort, color = cohort)) +
geom_point() +
geom_line() +
xlab("AgeInSamples")
ggplot(sParamsAll_3 |> filter(varNoIndex == "grIntT", mean > -5, mean < 5), aes(varNumeric1, mean, group = cohort, color = cohort)) +
geom_point() +
geom_line() +
xlab("AgeInSamples")
ggplot(sParamsAll_3 |> filter(varNoIndex == "phiBetaT", mean > -5, mean < 5), aes(varNumeric2, mean, group = cohort, color = cohort)) +
geom_point() +
geom_line() +
xlab("AgeInSamples") +
facet_wrap((~varNumeric1))
ggplot(sParamsAll_3 |> filter(varNoIndex == "sumWeightZ01"), aes(varNumeric1, mean, group = cohort, color = cohort)) +
geom_point() +
geom_line() +
xlab("AgeInSamples")
ggplot(sParamsAll_3 |> filter(varNoIndex == "sumZ01"), aes(varNumeric1, mean, group = cohort, color = cohort)) +
geom_point() +
geom_line() +
xlab("AgeInSamples")
Need to get rid of squared term for cohort-specific models
phiInt_3 <- sParamsAll_3 |> filter(varNoIndex == "phiIntTOut") |>
dplyr::select(mean, varNumeric1, cohort) |>
rename(phiIntOut = mean, ais = varNumeric1)
phiBetas_3 <- sParamsAll_3 |> filter(varNoIndex == "phiBetaT", mean > -5, mean < 5) |>
dplyr::select(mean, varNumeric1, varNumeric2, cohort) |>
pivot_wider(names_from = varNumeric1, values_from = mean, names_prefix = "beta") |>
rename(ais = varNumeric2)
phiSize_3 = expand.grid(weight = seq(-2,2,0.25), ais = 1:11, cohort = cohorts) |>
left_join(phiInt_3) |>
left_join(phiBetas_3) |>
mutate(pSurv = ilogit(phiIntOut + beta1*weight + beta2*weight^2))
ggplot(phiSize_3, aes(weight, pSurv, color = factor(cohort), group = cohort)) +
geom_line() +
geom_point() +
facet_wrap(~ais)